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1.  Introduction.  In  this  paper,  we  study  the  Local  Discontinuous  Galerkin 
(LDG)  methods  for  nonlinear,  convection-diffusion  systems  of  the  form 

dtVi  +  V  ■  F(u,  jD  u)  =  0,  in  (0,  T)  x  fi, 

where  L!  C  and  u  =  (ui, . . . ,  Um)*-  The  LDG  methods  are  an  extension  of  the 
Runge-Kutta  Discontinuous  Galerkin  (RKDG)  methods  for  the  nonlinear  hyper¬ 
bolic  system 

dtu  -f-  V  •  f  (u)  =  0,  in  (0,  T)  x  Q, 

introduced  by  the  authors  [12,13,14,15,16],  and  further  developed  by  Atkins  and 
Shu  [2],  Bassi  and  Rebay  [4],  Bey  and  Oden  [7],  Biswas,  Devine,  and  Flaherty  [8], 
deCougny  et  al.  [17],  Devine  et  al.  [19,  20],  Lowrie,  Roe  and  van  Leer  [30],  and  by 
Ozturan  et  al.  [33].  The  RKDG  methods  are  constructed  by  applying  the  explicit 
time  discretizations  introduced  by  Shu  [37]  and  Shu  and  Osher  [38,39]  to  a  space 
discretization  that  uses  discontinuous  basis  functions.  Since  the  space  discretization 
is  highly  local  in  character  and  produces  easily  invertible,  block-diagonal  mass 
matrices  and  since  the  time-marching  scheme  is  explicit,  the  RKDG  method  is  a 
highly  parallelizable  method;  see  Biswas,  Devine,  and  Flaherty  [8].  Moreover,  it  is 
not  only  a  formally  high-order  accurate  method  that  can  easily  handle  complicated 
geometries,  but  it  satisfies  a  ceU  entropy  inequality  that  enforces  a  nonlinear  L^- 
stability  property  even  without  the  slope  limiters  typical  of  this  method;  see  Jiang 
and  Shu  [27]. 

Extensions  of  the  RKDG  method  to  hydrodynamic  models  for  semiconductor 
device  simulation  have  been  constructed  by  Chen  et  al.  [9] ,  [10] .  In  these  extensions, 
approximations  of  the  derivatives  of  the  discontinuous  approximate  solution  are 
obtained  by  using  a  simple  projection  into  suitable  finite  elements  spaces.  This 
projection  requires  the  inversion  of  global  mass  matrices,  which  in  [9]  and  [10]  are 
‘lumped’  in  order  to  maintain  the  high  parallelizability  of  the  method.  Since  in  [9] 
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and  [10]  polynomials  of  degree  one  are  used,  the  ‘mass  lumping’  is  justified;  however, 
if  polynomials  of  higher  degree  were  used,  the  ‘mass  lumping’  needed  to  enforce 
the  full  parallelizability  of  the  method  could  cause  a  degradation  of  the  formal 
order  of  accuracy.  Fortunately,  this  is  not  an  issue  with  the  methods  proposed  by 
Bassi  and  Rebay  [5]  (see  also  Bassi  et  al  [6])  for  the  compressible  Navier-Stokes 
equations.  In  these  methods,  the  original  idea  of  the  RKDG  method  is  applied 
to  both  u  and  Du  which  are  now  considered  as  independent  unknowns.  Like  the 
RKDG  methods,  the  resulting  methods  are  highly  parallelizable  methods  of  high- 
order  accuracy  which  are  very  efficient  for  time-dependent,  convection-dominated 
flows.  The  LDG  methods  are  a  generalization  of  these  methods. 

The  basic  idea  to  construct  the  LDG  methods  is  to  suitably  rewritethe  convection- 
diffusion  system  into  a  larger,  degenerate,  first-order  system  and  then  discretize  it 
by  the  RKDG  method.  By  a  careful  choice  of  this  rewriting,  nonlinear  stability  can 
be  achieved  even  without  slope  limiters,  just  as  the  RKDG  method  in  the  purely 
hyperbolic  case;  see  Jiang  and  Shu  [27].  In  the  linear  case,  the  stability  result  leads 
to  the  sub-optimal  rate  of  (Arc)*’  for  the  L°°(0,  T;L2)-norm  of  the  error  if  polynomi¬ 
als  of  degree  at  most  k  are  used.  However,  these  estimates  are  sharp,  as  numerical 
evidence  reported  in  Bassi  et  al.  [6]  and  in  this  paper  indicate.  In  the  purely 
hyperbolic  case,  the  rate  of  convergence  of  is  recovered,  as  expected. 

Indeed,  this  is  the  same  rate  of  convergence  obtained  for  the  original  Discontinuous 
Galerkin  method  (introduced  by  Reed  and  Hill  [35])  for  purely  hyperbolic  case  by 
Johnson  and  Pitkaranta  [28]  and  confirmed  to  be  optimal  by  Peterson  [34],  LeSaint 
and  Raviait  [29]  proved  a  rate  of  convergence  of  (Acc)^  for  general  triangulations 
and  of  (Aa;)^'''^  for  Cartesian  grids;  Richter  [36]  obtained  the  optimal  rate  of  con¬ 
vergence  of  (Aa:)*'^^  for  some  structured  two-dimensional  non-Cartesian  grids.  The 
technique  for  proving  the  error  estimates  used  in  this  paper  is  different  from  the 
techniques  used  in  the  above  mentioned  papers.  It  is  very  simple  and  relies,  as  ex- 
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pected,  on  a  straightforward  combination  of  (i)  the  L^-stability  of  the  LDG  method 
and  of  (ii)  the  approjdmation  properties  of  the  finite  element  spaces. 

The  LDG  methods  introduced  in  this  paper  are  very  different  from  the  so-called 
Discontinuous  Galerkin  (DG)  method  for  parabolic  problems  introduced  by  Jamet 
[26]  and  studied  by  Eriksson,  Johnson,  and  Thom4e  [25],  Eriksson  and  Johnson  [21, 
22,  23,  24],  and  more  recently  by  Makridakis  and  Babuska  [31].  In  the  DG  method, 
the  approximate  solution  is  discontinuous  only  in  time,  not  in  space;  in  fact,  the 
space  discretization  is  the  standard  Galerkin  discretization  with  continuous  finite 
elements.  This  is  in  strong  contrast  with  the  space  discretizations  of  the  LDG 
methods  which  use  discontinuous  finite  elements.  To  emphasize  this  difference,  we 
call  the  methods  developed  in  this  paper  the  Local  Discontinuous  Galerkin  meth¬ 
ods.  We  also  must  emphasize  that  the  large  number  of  degrees  of  freedom  and  the 
restrictive  conditions  of  the  size  of  the  time  step  for  explicit  time-discretizations, 
render  the  LDG  methods  inefficient  for  diffusion-dominated  problems;  in  this  sit¬ 
uation,  the  use  of  methods  with  continuous-in-space  approximate  solutions  is  rec¬ 
ommended.  However,  as  for  the  successful  RKDG  methods  for  purely  hyperbolic 
problems,  the  extremely  local  domain  of  dependency  of  the  LDG  methods  allows 
a  very  efiicient  parallelization  that  by  far  compensates  for  the  extra  number  of 
degrees  of  freedom  in  the  case  of  convection-dominated  flows. 

Many  researchers  have  worked  in  the  devising  and  analysis  of  numerical  meth¬ 
ods  for  convection-dominated  flows.  In  particular,  Dawson  [18]  and,  more  recently, 
Arbogast  and  Wheeler  [1]  have  developed  and  analyzed  methods  that  share  several 
properties  with  the  LDG  methods:  They  use  discontinuous-in-space  approxima¬ 
tions,  they  are  locally  conservative,  and  they  approximate  the  diffusive  fluxes  with 
independent  variables  (by  using  a  mixed  method).  We  refer  the  reader  interested  in 
numerical  methods  for  convection-dominated  flows  to  [18]  and  [1]  and  the  references 


therein. 
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Another  numerical  method  that  uses  discontinuous  approximations  is  the  one 
proposed  and  studied  by  Baker  et  al.  [3].  This  method,  however,  is  not  for 
convection-dominated  flows  but  for  the  Stokes  problem.  The  advantage  of  using 
discontinuous  approximations  in  this  setting  is  that  this  allows  for  a  pointwise  ver¬ 
ification  of  the  incompressibility  condition  at  the  interior  of  each  triangle.  Optimal 
estimates  are  obtained. 

In  this  paper,  we  restrict  ourselves  to  the  semidiscrete  LDG  methods  for  convection- 
diffusion  problems  with  periodic  boundary  conditions.  Our  aim  is  to  clearly  display 
the  most  distinctive  features  of  the  LDG  methods  in  a  setting  as  simple  as  possible. 
The  fully  discrete  methods  for  convection-diffusion  problems  in  bounded  domains 
will  be  treated  in  a  forthcoming  paper.  This  paper  is  organized  as  follows;  In  §2, 
we  introduce  the  LDG  methods  for  the  simple  one-dimensional  case  d  =  1  in  which 

F(«,  Du)  =  f{u)  —  a{u)  dxU, 

u  is  a  scalar  and  a{u)  >  0  and  show  some  preliminary  numerical  results  displaying 
the  performance  of  the  method.  In  this  simple  setting,  the  main  ideas  of  how  to 
devise  the  method  and  how  to  analyze  it  can  be  clearly  displayed  in  a  simple  way. 
Thus,  the  L^-stability  of  the  method  is  proven  in  the  general  nonlinear  case  and 
the  rate  of  convergence  of  (Arc)^  in  the  L°°(0,  T;L^)-norm  for  polynomials  of  degree 
A:  >  0  in  the  linear  case  is  obtained;  this  estimate  is  sharp.  In  §3,  we  extend  these 
results  to  the  case  in  which  u  is  a  scalar  and 

Fi{u,Du)  =  fi{u)  -  ^  aij{u)dxjU, 
i<j<d 

where  defines  a  positive  semidefinite  matrix.  Again,  the  L^-stability  of  the 
method  is  proven  for  the  general  nonlinear  case  and  the  rate  of  convergence  of 
(Ax)''  in  the  L°°(0,r;L^)-norm  for  polynomials  of  degree  fc  >  0  and  arbitrary  tri¬ 
angulations  is  proven  in  the  linear  case.  In  this  case,  the  multidimensionality  of  the 
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problem  and  the  arbitrariness  of  the  grids  increase  the  technicality  of  the  analysis 
of  the  method  which,  nevertheless,  uses  the  same  ideas  of  the  one-dimensional  case. 
In  §4,  the  extension  of  the  LDG  method  to  multidimensional  systems  is  briefly 
described  and  concluding  remarks  are  offered. 

2.  The  LDG  methods  for  the  one-dimensional  case.  In  this  section,  we 
present  and  analyze  the  LDG  methods  for  the  following  simple  model  problem: 


dtu  +  dx  (/(«)  -  a(u)  dxu)  =  0 

in  (0,r)x  (0,1), 

(2.1a) 

II 

0 

II 

0 

on  (0, 1), 

(2.1b) 

with  periodic  boundary  conditions. 

a.  General  formulation  and  main  properties.  To  define  the  LDG  method. 

we  introduce  the  new  variable  q  =  s/  o.{u)  dx  u 

and  rewrite  the  problem  (2.1)  as 

follows: 

dtu  +  dx  {f{u)  -  -sj a{u)  q)  =  0 

in  (0,T)x(0,l), 

(2.2a) 

q-dx  g{u)  =  0 

in  (0,T)x(0,l), 

(2.2b) 

u{t  =  0)  =  wo) 

on  (0,1), 

(2.2c) 

where  g{u)  =  /“  a{s)  ds.  The  LDG  method  for  (2.1)  is  now  obtained  by  simply 
discretizing  the  above  system  with  the  Discontinuous  Galerkin  method. 

To  do  that,  we  follow  [13]  and  [14].  We  define  the  flux  h  =  ( hq  )*  as  follows: 

h(u,  q)  =  ( f{u)  -  a/ a{u)  q ,  -g{u)  )*.  (2.3) 

For  each  partition  of  the  interval  (0, 1),  {xj+i/2}^0’  h  ~  0+1/2) ■> 

and  Aajj  =  for  j  =  1, . . . ,  N]  we  denote  the  quantity  maxi<j<iv 

by  Ar  .  We  seek  an  approximation  =  {uh^qh)*  to  w  =  {u,q)*  such  that  for 
each  time  t  E  [0,T],  both  Uh{t)  and  qh{t)  belong  to  the  finite  dimensional  space 

14  =  =  {^  G  L\0, 1)  :  v\r,  G  j  =  l,...,N},  (2.4) 
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where  P''(7)  denotes  the  space  of  polynomials  in  I  of  degree  at  most  k.  In  order 
to  determine  the  approximate  solution  {uh,  qh)>  we  first  note  that  by  multiplying 
(2.2a),  (2.2b),  and  (2.2c)  by  arbitrary,  smooth  functions  Vu,  Vg,  anduj,  respectively, 
and  integrating  over  Ij,  we  get,  after  a  simple  formal  integration  by  parts  in  (2.2a) 
and  (2.2b), 

/  dtu{x,t)vu{x)dx-  /  hu{w{x,t))dxVuix)dx 

Jh 

+  K{^{Xj+i/2,t))Vu{x-T^y2)  -  hn{^{Xj_i/2,t))Vu{x-f_y^)  =  0, 

(2.5a) 

/  q{x,t)vq{x)dx  -  /  hg{w{x,t))dxVq{x)dx 

Jij  Jij 

+  hq{^N{Xj+il2,t))Vq{x-,^)  -  hq{^{x^-il2,t))  Vq{x^_^f^)  =  0, 

(2.5b) 

/  u{x,Q)vi{x)dx  —  /  uo{x)vi{x)  dx.  (2.5c) 

Next,  we  replace  the  smooth  functions  Vu,  Vq,  and  Vi  by  test  functions  Vh,u,  Vh,q,  and 
Vh,i\  respectively,  in  the  finite  element  space  Vh  and  the  exact  solution  w  =  {u,  qY 
by  the  approximate  solution  Wh  =  {uh,qhY-  Since  this  function  is  discontinuous 
in  each  of  its  components,  we  must  also  replace  the  nonlinear  fiux  h(w(a::j+i/2,f)) 
by  a  numerical  flux  h(w)j+i/2(i)  =  {hu{'Wh)j+i/2{t),hq(wh)j+i/2{i))  that  will  be 
suitably  chosen  later.  Thus,  the  approximate  solution  given  by  the  LDG  method 
is  defined  as  the  solution  of  the  following  weaJc  formulation: 
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It  only  remains  to  choose  the  numerical  flux  h{wh}j+i/2{'t)-  We  use  the  notation: 
[p]=j)+-p-,  and  p=^ip'^  +P~), 

and  pf+if2  =  P(®*+i/2)‘  consistent  with  the  type  of  numerical  fluxes  used  in 
the  RKDG  methods,  we  consider  numerical  fluxes  of  the  form 

h{-Wh)j+i/2{t)  =  H^h{xJ+i/2,'t),^h{x^+i/2>i))> 

that  (i)  are  locally  Lipschitz  and  consistent  with  the  flux  h,  (ii)  allow  for  a  local 
resolution  of  qh  in  terms  of  Uh,  (iii)  reduce  to  an  E-flux  (see  Osher  [32])  when 
a(-)  =  0,  and  that  (iv)  enforce  the  L^-stability  of  the  method. 

To  reflect  the  convection-diffusion  nature  of  the  problem  under  consideration, 
we  write  our  numerical  flux  as  the  sum  of  a  convective  flux  and  a  diffusive  flux: 

h(w“,w+)  =  hcon«(w“,W'‘)-|-h^ji//(w“,w+).  (2.7a) 


The  convective  flux  is  given  by 

hcon«(w",w+)  =  (/(u“,  «■*■),  0)*,  (2.7b) 

where  f{u~,u^)  is  any  locally  Lipschitz  E-flux  consistent  with  the  nonlinearity  /, 
and  the  diffusive  flux  is  given  by 

hdi//(w~,w+)  =  (  -  -9{u)  y  -  ^diff  [w],  (2.7c) 

where 


T)> 

(2.7d) 

Ci2  =  Ci2(w“,  w"*")  is  locally  Lipschitz, 

(2.7e) 

ci2  =  0  when  a(’)  =  0. 

(2.7f) 

We  claim  that  this  flux  satisfies  the  properties  (i)  to  (iv). 
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Let  us  prove  our  claim.  That  the  flux  h  is  consistent  with  the  flux  h  easily 
follows  from  their  definitions,  (2.3)  and  (2.7).  That  h  is  locally  Lipschitz  follows 
from  the  fact  that  /(•,  •)  is  locally  Lipschitz  and  from  (2.7d);  we  assume  that  /(•) 
and  o(-)  are  locally  Lipschitz  functions,  of  course.  Property  (i)  is  hence  satisfied. 

That  the  approximate  solution  qn  can  be  resolved  element  by  element  in  terms  of 
Uh  by  using  (2.6b)  follows  from  the  fact  that,  by  (2.7c),  the  flux  hq  =  -g{u)-ci2  [u] 
is  independent  of  qn.  Property  (ii)  is  hence  satisfied. 

Property  (iii)  is  also  satisfied  by  (2.7f)  and  by  the  construction  of  the  convective 
flux. 

To  see  that  the  property  (iv)  is  satisfied,  let  us  first  rewrite  the  flux  h  in  the 
following  way: 


h(w  ,  w'*')  =  ( 


IMI 

[u] 


-g{u)  )*-C[w], 


where 


with  0(u)  defined  by  0(u)  =  /“  f{s)  ds.  Since  /(•,  •)  is  an  E-flux, 

^11  =  7^2  /  {f{s)- f{u~,u-^))ds>0, 

Ju> 

and  so,  by  (2.7d),  the  matrix  C  is  semipositive  definite.  The  property  (iv)  follows 
from  this  fact  and  from  the  following  result. 


Proposition  2.1.  (L^-stability)  We  have, 


1 

2 


J  ul{x,T)dx  +  ql{x,t)dxdt+eT,c{[v^h])  ul{x)dx, 


©T,c([w/i])  =  /  J]  I  [Wh(t)]*c  [whit)]  I  dt. 

"'O  l<j<N  ^  ^  i+1/2 


where 
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This  result  will  be  proven  in  section  §2.c.  Thus,  this  shows  that  the  flux  h  given 
by  (2.7)  does  satisfy  the  properties  (i)  to  (iv)-  as  claimed. 

Now,  we  turn  to  the  question  of  the  quality  of  the  approximate  solution  defined 
by  the  LDG  method.  In  the  linear  case  /'  =  c  and  a(-)  =  a,  from  the  above  stability 
result  and  from  the  the  approximation  properties  of  the  finite  element  space  14, 
we  can  prove  the  following  error  estimate.  We  denote  the  L^(0,  l)-norm  of  the  ^-th 
derivative  of  «  by  |  tt  fy 

Theorem  2.2.  (L^-error  estimate)  Let  e  be  the  approximation  error  w—wh-  Then 
we  have, 

\eu{x,T)\'^dx  + \eg{x,  t)  pda;dt  +  ©T,c([e])  I  <  C{Ax)^, 

where  C  =  C{k,  |  u  4+1,  |  n  4+2)-  In  the  purely  hyperbolic  case  a  =  0,  the  constant 
C  is  of  order  (Ax)^/^  and  in  the  purely  parabolic  case  c  =  0,  the  constant  C  is  of 
order  Ax  for  even  values  of  k  for  uniform  grids  and  for  C  identically  zero. 

This  result  will  be  proven  in  section  §2.d.  The  above  error  estimate  gives  a 
suboptimal  order  of  convergence,  but  it  is  sharp  for  the  LDG  methods.  Indeed, 
Bassi  et  al  [6]  report  an  order  of  convergence  of  order  k  +  1  for  even  values  of 
k  and  of  order  k  for  odd  values  of  k  for  a  steady  state,  purely  elliptic  problem 
for  uniform  grids  and  for  C  identically  zero.  Our  numerical  results  for  a  purely 
parabolic  problem  give  the  same  conclusions;  see  Table  5  in  the  section  §2.b. 

Our  error  estimate  is  also  sharp  in  that  the  optimal  order  of  convergence  of 
fc  +  1/2  is  recovered  in  the  purely  hyperbolic  case,  as  expected.  This  improvement 
of  the  order  of  convergence  is  a  reflection  of  the  semipositive  definiteness  of  the 
matrix  C,  which  enhances  the  stability  properties  of  the  LDG  method.  Indeed, 
since  in  the  purely  hyperbolic  case 

©T,c([w^,])=  /  {  [«a(t)]*cii  K(t)]  1  dt, 

*^0  l<j<N  ^  Jj+1/2 
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the  method  enforces  a  control  of  the  jumps  of  the  variable  Uh,  as  shown  in  Propo¬ 
sition  2.1.  This  additional  control  is  reflected  in  the  improvement  of  the  order  of 
accuracy  from  k  in  the  general  case  to  A;  -|-  1/2  in  the  purely  hyperbolic  case. 

However,  this  can  only  happen  in  the  purely  hyperbolic  case  for  the  LDG  meth¬ 
ods.  Indeed,  since  Cn  =  0  for  c  =  0,  the  control  of  the  jumps  of  Uh  is  not  enforced 
in  the  purely  parabolic  case.  As  indicated  by  the  ninnerical  experiments  of  Bassi  et 
al.  [6]  and  those  of  section  §2.b  below,  this  can  result  in  the  effective  degradation  of 
the  order  of  convergence.  To  remedy  this  situation,  the  control  of  the  jumps  of  Uh 
in  the  purely  parabolic  case  can  be  easily  enforced  by  letting  Cn  be  strictly  positive 
if  |c|  -h  |a|  >  0.  Unfortunately,  this  is  not  enough  to  guarantee  an  improvement 
of  the  accuracy:  an  additional  control  on  the  jumps  of  Qh  is  required!  This  can  be 
easily  achieved  by  allowing  the  matrix  C  to  be  symTnetric  arid  positive  definite  when 
o  >  0.  In  this  case,  the  order  of  convergence  of  fc  -t- 1/2  can  be  easily  obtained  for 
the  general  convection-diffusion  case.  However,  this  would  force  the  matrix  entry 
C22  to  be  nonzero  and  the  property  (ii)  of  local  resolvability  of  qn  in  terms  of  Uh 
would  not  be  satisfied  anymore.  As  a  consequence,  the  high  parallelizability  of  the 
LDG  would  be  lost. 

The  above  result  shows  how  strongly  the  order  of  convergence  of  the  LDG  meth¬ 
ods  depend  on  the  choice  of  the  matrix  C.  In  fact,  the  numerical  results  of  section 
§2.b  in  uniform  grids  indicate  that  with  yet  another  choice  of  the  matrix  C,  see 
(2.9),  the  LDG  method  converges  with  the  optimal  order  of  fc  -f  1  in  the  general 
case.  The  analysis  of  this  phenomenon  constitutes  the  subject  of  ongoing  work. 


b.  Preliminary  numerical  results. 

In  this  section  we  provide  preliminary  numerical  results  for  the  schemes  discussed 
in  this  paper.  We  will  only  provide  results  for  the  following  one  dimensional,  linear 
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convection  diffusion  equation 

dtu  +  cdxU  —  ad^u  =  0  in  (0, T)  x  (0, 2 tt), 

M(t  =  0,a;)  =  sin(a;),  on  (0,27r), 

where  c  and  a  >  0  are  both  constants;  periodic  boundary  conditions  are  used.  The 
exact  solution  is  u{t,x)  =  e““*sin(a:  —  ct).  We  compute  the  solution  up  to  T  =  2, 
and  use  the  LDG  method  with  C  defined  by 


We  notice  that,  for  this  choice  of  fluxes,  the  approximation  to  the  convective  term 
cux  is  the  standard  upwinding,  and  that  the  approximation  to  the  diffusion  term 
ad^uis  the  standard  three  point  central  difference,  for  the  case.  On  the  other 
hand,  if  one  uses  a  central  flux  corresponding  to  Ci2  =  -C21  =  0,  one  gets  a  spread- 
out  five  point  central  difference  approximation  to  the  diffusion  term  adlu. 

The  LDG  methods  based  on  with  fc  =  1,2,3,4  are  tested.  Elements  with 
equal  size  are  used.  Time  discretization  is  by  the  third-order  accurate  TVD  Runge- 
Kutta  method  [38],  with  a  sufficiently  small  time  step  so  that  error  in  time  is 
negligible  comparing  Avith  spatial  errors.  We  list  the  Loo  errors  and  numerical 
orders  of  accuracy,  for  Uh^  as  well  as  for  its  derivatives  suitably  scaled  Ax'^d^Uh 
for  1  <  m  <  fc,  at  the  center  of  of  each  element.  This  gives  the  complete  description 
of  the  error  for  Uh  over  the  whole  domain,  as  Uh  in  each  element  is  a  polynomial 
of  degree  k.  We  also  list  the  Loo  errors  and  numerical  orders  of  accuracy  for  qn  at 
the  element  center. 

In  aH  the  convection-diffusion  runs  with  a  >  0,  accuracy  of  at  least  (fc  -f  l)-th 
order  is  obtained,  for  both  Uh  and  9/1,  when  P^  elements  are  used.  See  Tables  1  to 
3.  The  P'*  case  for  the  purely  convection  equation  a  =  0  seems  to  be  not  in  the  as¬ 
ymptotic  regime  yet  with  iV  =  40  elements  (further  refinement  with  Ai  =  80  suffers 
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from  round-oflf  effects  due  to  our  choice  of  non-orthogonal  basis  functions),  Table 
4.  However,  the  absolute  values  of  the  errors  are  comparable  with  the  convection 
dominated  case  in  Table  3. 


Table  1.  The  heat  equation  a  =  1,  c  =  0.  Loo  errors  and  numerical  order  of 
accuracy,  measured  at  the  center  of  each  element,  for  Ax'^d^Uh  for  0  <  m  <  A:, 
and  for  qh- 


k 

variable 

N  =  10 

error 

N  = 

error 

20 

order 

AT  =  40 

error  order 

u 

4.55E-4 

5.79E-5 

2.97 

7.27E-6 

2.99 

1 

AxdxU 

9.01E-3 

2.22E-3 

2.02 

5.56E-4 

2.00 

Q 

4.17E-5 

2.48E-6 

4.07 

1.53E-7 

4.02 

u 

1.43E-4 

1.76E-5 

3.02 

2.19E-6 

3.01 

2 

AxdxU 

7.87E-4 

1.03E-4 

2.93 

1.31E-5 

2.98 

(Aa;)^  dxU 

1.64E-3 

2.09E-4 

2.98 

2.62E-5 

2.99 

q 

1.42E-4 

1.76E-5 

3.01 

2.19E-6 

3.01 

u 

1.54E-5 

9.66E-7 

4.00 

6.11E-8 

3.98 

AxdxU 

3.77E-5 

2.36E-6 

3.99 

1.47E-7 

4.00 

3 

(Aa;)2  dlu 

1.90E-4 

1.17E-5 

4.02 

7.34E-7 

3.99 

(Ax)^  dxU 

2.51E-4 

1.56E-5 

4.00 

9.80E-7 

4.00 

q 

1.48E-5 

9.66E-7 

3.93 

6.11E-8 

3.98 

u 

2.02E-7 

5.51E-9 

5.20 

1.63E-10 

5.07 

AxdxU 

1.65E-6 

5.14E-8 

5.00 

1.61E-9 

5.00 

4 

(Aa:)2  d^u 

6.34E-6 

2.04E-7 

4.96 

6.40E-9 

4.99 

{Axfd^u 

2.92E-5 

9.47E-7 

4.95 

2.99E-8 

4.99 

(Ax)^  d^u 

3.03E-5 

9.55E-7 

4.98 

2.99E-8 

5.00 

q 

2.10E-7 

5.51E-9 

5.25 

1.63E-10 

5.07 
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Table  2.  The  convection  diifusion  equation  a  =  1,  c  =  1.  Loo  errors  and  numerical 
order  of  accuracy,  measured  at  the  center  of  each  element,  for  Ax^d^  Uh  for  0  < 
m<  k,  and  for  qh- 


k 

variable 

II 

O 

N  = 

20 

II 

o 

error 

error 

order 

error 

order 

u 

6.47E-4 

1.25E-4 

2.37 

1.59E-5 

2.97 

1 

Axd^u 

9.61E-3 

2.24E-3 

2.10 

5.56E-4 

2.01 

<1 

2.96E-3 

1.20E-4 

4.63 

1.47E-5 

3.02 

u 

1.42E-4 

1.76E-5 

3.02 

2.18E-6 

3.01 

2 

AxdxU 

7.93E-4 

1.04E-4 

2.93 

1.31E-5 

2.99 

(Ax)^  d^u 

1.61E-3 

2.09E-4 

2.94 

2.62E-5 

3.00 

Q 

1.26E-4 

1.63E-5 

2.94 

2.12E-6 

2.95 

u 

1.53E-5 

9.75E-7 

3.98 

6.12E-8 

3.99 

AxdxU 

3.84E-5 

2.34E-6 

4.04 

1.47E-7 

3.99 

3 

(Ax)^  d'^u 

1.89E-4 

1.18E-5 

4.00 

7.36E-7 

4.00 

(Ax)3  dlu 

2.52E-4 

1.56E-5 

4.01 

9.81E-7 

3.99 

Q 

1.57E-5 

9.93E-7 

3.98 

6.17E-8 

4.01 

u 

2.04E-7 

5.50E-9 

5.22 

1.64E-10 

5.07 

AxdxU 

1.68E-6 

5.19E-8 

5.01 

1.61E-9 

5.01 

4 

(C^fdlu 

6.36E-6 

2.05E-7 

4.96 

6.42E-8 

5.00 

(Ax)^  dl.u 

2.99E-5 

9.57E-7 

4.97 

2.99E-8 

5.00 

(Ax)^  dlu 

2.94E-5 

9.55E-7 

4.95 

3.00E-8 

4.99 

q 

1.96E-7 

5.35E-9 

5.19 

1.61E-10 

5.06 
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Table  3.  The  convection  dominated  convection  diffusion  equation  a  =  0.01,c  =  l. 
Loo  errors  and  numerical  order  of  accuracy,  measured  at  the  center  of  each  element, 
for  Ax"^d^Uh  for  0  <  m  <  fc,  and  for  qh- 


k 

variable 

N  =  10 

N  = 

20 

o 

II 

error 

error 

order 

error 

order 

u 

7.14E-3 

9.30E-4 

2.94 

1.17E-4 

2.98 

1 

AxdxU 

6.04E-2 

1.58E-2 

1.93 

4.02E-3 

1.98 

Q 

8.68E-4 

1.09E-4 

3.00 

1.31E-5 

3.05 

u 

9.59E-4 

1.25E-4 

2.94 

1.58E-5 

2.99 

2 

AxdxU 

5.88E-3 

7.55E-4 

2.96 

9.47E-5 

3.00 

{Axfd^u 

1.20E-2 

1.50E-3 

3.00 

1.90E-4 

2.98 

Q 

8.99E-5 

l.llE-5 

3.01 

l.lOE-6 

3.34 

u 

l.llE-4 

7.07E-6 

3.97 

4.43E-7 

4.00 

AxdxU 

2.52E-4 

1.71E-5 

3.88 

1.07E-6 

4.00 

3 

(Aa:)2  d^u 

1.37E-3 

8.54E-5 

4.00 

5.33E-6 

4.00 

(Axfd^u 

1.75E-3 

1.13E-4 

3.95 

7.11E-6 

3.99 

Q 

1.18E-5 

7.28E-7 

4.02 

4.75E-8 

3.94 

u 

1.85E-6 

4.02E-8 

5.53 

1.19E-9 

5.08 

AxdxU 

1.29E-5 

3.76E-7 

5.10 

1.16E-8 

5.01 

4 

(Ax)2  dlu 

5.19E-5 

1.48E-6 

5.13 

4.65E-8 

4.99 

{Axfdlu 

2.21E-4 

6.93E-6 

4.99 

2.17E-7 

5.00 

(Ax)^  d^u 

2.25E-4 

6.89E-6 

5.03 

2.17E-7 

4.99 

Q 

3.58E-7 

3.06E-9 

6.87 

5.05E-11 

5.92 
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Table  4.  The  convection  equation  a  =  0,  c  =  1.  Loo  errors  and  numerical  order  of 
accuracy,  measured  at  the  center  of  each  element,  for  Ax^d^  for  0  <  m  <  k. 


k 

variable 

JV  =  10 

o 

(M 

11 

II 

o 

error 

error 

order 

error 

order 

1 

u 

7.24E-3 

9.46E-4 

2.94 

1.20E-4 

2.98 

AxdxU 

6.09E-2 

1.60E-2 

1.92 

4.09E-3 

1.97 

u 

9.96E-4 

1.28E-4 

2.96 

1.61E-5 

2.99 

2 

AxdxU 

6.00E-3 

7.71E-4 

2.96 

9.67E-5 

3.00 

(Aa:)2  dlu 

1.23E-2 

1.54E-3 

3.00 

1.94E-4 

2.99 

u 

1.26E-4 

7.50E-6 

4.07 

4.54E-7 

4.05 

3 

AxdxU 

1.63E-4 

2.00E-5 

3.03 

1.07E-6 

4.21 

(Aa;)2  dlu 

1.52E-3 

9.03E-5 

4.07 

5.45E-6 

4.05 

(Ax)3  dlu 

1.35E-3 

1.24E-4 

3.45 

7.19E-6 

4.10 

u 

3.55E-6 

8.59E-8 

5.37 

3.28E-10 

8.03 

Ax  dxU 

1.89E-5 

1.27E-7 

7.22 

1.54E-8 

3.05 

4 

(Aa:)^  dlu 

8.49E-5 

2.28E-6 

5.22 

2.33E-8 

6.61 

(Ax)^  dlu 

2.36E-4 

5.77E-6 

5.36 

2.34E-7 

4.62 

(Ax)^dlu 

2.80E-4 

8.93E-6 

4.97 

1.70E-7 

5.72 

Finally,  to  show  that  the  order  of  accuracy  could  really  degenerate  to  k  for  P^, 
as  was  already  observed  in  [6],  we  rerun  the  heat  equation  case  a  =  1,  c  =  0  with 
the  central  flux 


This  time  we  can  see  that  the  global  order  of  accuracy  in  Lqo  is  only  k  when  is 
used  with  an  odd  value  of  k. 
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Table  5.  The  heat  equation  a  =  1,  c  =  0.  Loo  errors  and  numerical  order  of 
accuracy,  measured  at  the  center  of  each  element,  for  Ax'^d^Uh  for  0  <  m  <  fc, 


and  for  qh,  using  the  central  flux. 


k 

variable 

AT  =  10 

N  = 

20 

o 

11 

error 

error 

order 

error 

order 

u 

3.59E-3 

8.92E-4 

2.01 

2.25E-4 

1.98 

1 

Axd^u 

2.10E-2 

1.06E-2 

0.98 

5.31E-3 

1.00 

Q 

2.39E-3 

6.19E-4 

1.95 

1.56E-4 

1.99 

u 

6.91E-5 

4.12E-6 

4.07 

2.57E-7 

4.00 

2 

AxdxU 

7.66E-4 

1.03E-4 

2.90 

1.30E-5 

2.98 

2.98E-4 

1.68E-5 

4.15 

1.03E-6 

4.02 

Q 

6.52E-5 

4.11E-6 

3.99 

2.57E-7 

4.00 

u 

1.62E-5 

l.OlE-6 

4.00 

6.41E-8 

3.98 

AxdxU 

1.06E-4 

1.32E-5 

3.01 

1.64E-6 

3.00 

3 

iAx)^d^u 

1.99E-4 

1.22E-5 

4.03 

7.70E-7 

3.99 

iAxfd^u 

6.81E-4 

8.68E-5 

2.97 

1.09E-5 

2.99 

Q 

1.54E-5 

l.OlE-6 

3.93 

6.41E-8 

3.98 

u 

8.25E-8 

1.31E-9 

5.97 

2.11E-11 

5.96 

AxdxU 

1.62E-6 

5.12E-8 

4.98 

1.60E-9 

5.00 

4 

{Axf  dlu 

1.61E-6 

2.41E-8 

6.06 

3.78E-10 

6.00 

(Ax)^  d^u 

2.90E-5 

9.46E-7 

4.94 

2.99E-8 

4.99 

{Ax^d^u 

5.23E-6 

7.59E-8 

6.11 

1.18E-9 

6.01 

Q 

7.85E-8 

1.31E-9 

5.90 

2.11E-11 

5.96 

c.  Proof  of  the  nonlinear  stability.  In  this  section,  we  prove  the  the  nonlinear 
stability  result  of  Proposition  2.1.  To  do  that,  we  first  show  how  to  obtain  the 
corresponding  stability  result  for  the  exact  solution  and  then  mimic  the  argument 
to  obtain  Proposition  2.1. 

The  continuous  case  as  a  model.  We  start  by  rewriting  the  equations  (2.5a) 
and  (2.5b),  in  compact  form.  If  in  equations  (2.5a)  and  (2.5b)  we  replace  Vu{x)  and 
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analyzes  proposed  legislation  and  regulations  for  their  effects  on  the 
industry. 

•  Provides  actuarial  and  statistical  foundations  for  HIAA  policy. 

•  Develops  policy  proposals  based  on  members'  input  and  government 
and  industry  data. 

Public  Affairs 

•  Conducts  advocacy  campaigns  nationwide  and  in  the  states, 
coordinating  with  HIAA  loftiyists,  member  companies,  and  allies. 

•  Builds  (and  activates)  coalition  and  grassroots  and 
"grasstops"networks. 

•  Evaluates  political  climate;  develops  and  tests  messages  for  advocacy 
and  product  support;  uses  public  opinion  research  (surveys,  polling, 
and  focus  groups). 

•  Commimicates  with  media  and  public  through  press  conferences, 
articles,op-eds,  letters  to  the  editor,  advertising,  and  other  means. 

•  Conducts  media  tours,  places  speakers,  offers  media  training  to 
members. 

Information 

•  Provides  members  access  to  federal  and  state  legislative  and 
regulatory  information  and  to  HIAA  reports  and  information  through 
an  online  service  called  HiWire. 

•  Provides  information  to  consumers,  the  media,  insurance 
professionals  and  others  via  the  World  Wide  Web 
(http://www.hiaa.org). 

•  Maintains  a  library  of  health  insurance  related  material  of  interest  to 
researchers  and  industry  professionals. 

Publications 

•  HIAA  publishes  books,  white  papers,  reports,  studies,  promotional 
materials,  buyers'  guides,  lobbying  kits,  and  other  materials. 

•  The  Source  Book  of  Health  Insurance  Data,  our  flagship  publication, 
is  used  all  over  the  world  by  readers  interested  in  health  care  policy. 

•  The  HIAA  Reporter,  a  bimonthly  newsletter,  keeps  members  and 
interested  parties  abreast  of  industry  and  Association  activities. 

Insurance  Education 

•  Offers  innovative  coursework  by  mail  on  every  aspect  of  health 
insurance. 
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Health  Insurance  Association  of  America  -  Purpose  and  Programs 


http  ://www  .hiaa.org/abouthiaa/purpose.html 


•  Enrolls  over  30,000  students  a  year. 

•  Offers  professional  designations  and  CEU  credits. 

Prevailing  Healthcare  Charges  System  (PHCS) 

•  Preeminent  source  of  data  on  health  care  charges,  based  on  millions 
of  claims  records  and  used  by  hundreds  of  insurers. 

•  Offers  data  by  subscription,  in  all  formats  and  online. 

Meetings 

•  In  addition  to  on-going  Board  and  committee  meetings,  HIAA  holds 
seminars,  workshops,  and  fly-ins  on  special  subjects. 

•  HIAA  holds  an  annual  conference,  the  Insurance  Forum,  which  is 
open  to  members  and  non-members  and  features  an  exhibit  hall  and 
sessions  on  public  policy,  products,  markets,  and  other  industry 
issues. 


About  HIAA  Links 

Purpose  &  Programs  |  Vision  for  Healthcare  Reform 
Code  of  Ethics  |  Membership  Information 
HIAA’s  Members 


Whats  New  I  About  HIAA  I  Newsroom  I  Consumer  Information 
HIAA  Publications  [  Insurance  Education  Program 
Prevailing  Healthcare  Charges  System 
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Vq{x)  by  Vuix.t)  and  Vq{x,t),  respectively,  add  the  resulting  equations,  sum  on  j 
from  1  to  N ^  and  integrate  in  time  from  0  to  T,  we  obtain  that 

5(w,  v)  =  0,  V  smooth  v(t),  Vt  G  (0,T),  (2.10a) 


where 


B{w,v)=  [  [  dtu{x,t)vu{x,t)dxdt+  [  [  q{x,t)vq{x,t)  dx  dt 

Jo  Jo  Jo  Jo 

-  [  [  \i{-w{x^t)Y  ^x^^{x,t)  dxdt^  (2.10b) 

Jo  Jo 

using  the  fact  that  h(w(a;,  t)y  dx^{x,  t)  =  dx{(f>{u)  -  g{u)  g )  is  a  complete  deriva¬ 
tive,  we  see  that 

B{w,w)  =  -  /  u^{x,T)dx+  /  /  q^{x,t)dxdt-  -  /  ul{x)dx, 

2  Jo  Jo  Jo  2  Jo  (2.11) 

and  that  H(w,  w)  =  0,  by  (2.10a),  we  immediately  obtain  the  following  L^-stability 

result: 


This  is  the  argument  we  have  to  mimic  in  order  to  prove  Proposition  2.1. 

The  discrete  case.  Thus,  we  start  by  finding  a  compact  form  of  equations 
(2.6a)  and  (2.6b).  If  we  replace  Vh,u{x)  and  Vh,qix:)  by  Vh,u{x,t)  and  Vh,q{x,t)  in 
the  equations  (2.6a)  and  (2.6b),  add  them  up,  sum  on  j  from  1  to  iV  and  integrate 
in  time  from  0  to  T,  we  obtain 


-B;.(w;„v;,)  =  0,  yvh{t)eV^xV^,  VfG(0,r).  (2.12a) 


where 


( w/j, ,  )  —  [  [  dx  dt  V}i^q{^x^f)  dx  dt 

Jo  Jo  do  Jo 

-[  H^hYj+i/2{'t)[Mi)]j+i/2dt 

■^0  l<j<N 

-  j  Yj  [  ^{^h{x,t)y  dxVh{x,t)dxdt.  (2.12b) 


Next,  we  obtain  an  expression  for  Bh{'Wh,'Wh)-  It  is  contained  in  the  following 


result. 


Lemma  2.3.  l^e  have 


Bh{wh,Wh)  ul{x,T)dx  +  ^{x,t)dxdt+eT,c{[wh])-^  ul{x,0)dx, 

where  ©r,c([wh])  is  defined  in  Proposition  2.1. 

Next,  since  Bh{y^h,'^h)  =  0,  by  (2.12a),  we  get  the  inequality 

-  ul{x,T)dx  +  ql{x,t)dxdt  +  eT,c{[^h])^  ul{x,0)dx 

from  which  Proposition  2.1  easily  follows,  since 

ul{x,0)dx  ul{x)dx, 

by  (2.5c).  It  remains  to  prove  Lemma  2.3. 

Proof  Lemma  2.3.  After  setting  Vh  =  Wft,  in  (2.12b),  we  get 

B{wh,Wh)  [  ul{x,T)dx+  f  f  ql{x,t)dxdt+  f  &diss{t)  dt  -  ^  [  ul{x,0)dx, 

^  Jo  Jq  Jo  Jo  ^  Jo 


where 

®dissit)  =  -  ^  |h(w/,)5^i/2(t)  [’Wh(t)]j+i/2+  /  h(wh{x,t)YdxWh{x,t)dxY 

i<j<N  ^  J 

It  only  remains  to  show  that  &diss{t)  dt  =  ©T,c([wh,]). 

To  do  that,  we  proceed  as  follows.  Since 

h(w;i(a:,  t)Y  dx  Wh {x,  t)  =  ( f{uh)  -  \/a{uh)  qn  )dxUh-  g{uh)  dx  qn 

r'^h 

=  9x{  f{s)  ds  -  g{uh)  qh) 

=  dx  {d>{uh)  -  g{uh)  qh) 


=  dxH{wh{x,t)), 
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we  get 


^dissi^^  E  {  [H{-Wh{t))]j+i/2  -  h(Wh)5+i/2W  [Wh(i)]j+l/2| 
l<j<N  ^  ^ 

-  Y,  {  “Mw^)*(t)[w;,(t)]| 

1  <ro<r  AT  V  J  j~ 


l<j<N 

Since,  by  the  definition  of  H, 


'  J+l/2 


[H{wh{t))]  =  [  ^{Uh{t))]-  [giuhit))  qh{t)] 

=  [^(«/i(t))]  -  [9{uh{t))]qh{t)  -  [Qh{t)]g{uh{t)), 


and  since  {hu,  hqf  —  h,  we  get 

®di8s{t)=  V  \[(f>{uh{t))]-[g{uh{t))]qh{t)-[uh{t)]hu\ 

l<j<N  ^  Jj+l/S 

+  Y  (  “  [9/i(i)]5(«h)(*)- ^9}  » 

This  is  the  crucial  step  to  obtain  the  l!^ -stability  of  the  LDG  methods,  since  the 
above  expression  gives  us  key  information  about  the  form  that  the  flux  h  should 
have  in  order  to  make  &diss{i)  a  nonnegative  quantity  and  hence  enforce  the  bi¬ 
stability  of  the  LDG  methods.  Thus,  by  taking  h  as  in  (2.7a),  we  get 

&dissit)  =  (  [wh(t)]*C  [Whit)]  \ 

i<i<^  ^  ^  i+1/2 

and  the  result  follows.  This  completes  the  proof.  □ 


This  completes  the  proof  of  Proposition  2.1. 

c.  The  error  estimate  in  the  linear  case.  In  this  section,  we  prove  the  error 
estimate  of  Theorem  2.2  which  holds  for  the  linear  case  /'(•)  =  c  and  a(-)  =  a.  To  do 
that,  we  fibrst  show  how  to  estimate  the  error  between  the  solutions  w^,  =  {ui,,  qi,)*, 
i'  =  1, 2,  of 


dt  Ui,  +  dx  (/(«,/)  -  \/a('u„)  qv)  =0  in  (0,  T)  x  (0, 1), 
qv  -  dxgiu^,)  =  0 
U[,{t  =  0)  = 


in(0,T)x(0,l), 
on  (0, 1). 


20 


Then,  we  mimic  the  argument  in  order  to  prove  Theorem  2.2. 

The  continuous  case  as  amodel.  By  the  definition  of  the  form  (•,■),  (2.10b), 
we  have,  for  i/  =  1, 2, 

B{wv,y)  =  0,  V  smooth  v(t),  Vt  €  (0,T). 

Since  in  this  case,  the  form  is  bilinear,  from  the  above  equation  we  obtain 

the  so-called  error  equation: 

H(e,v)  =  0,  V  smooth  v(t),  Vt  G  (0,T), 


where  e  =  Wi  —  W2.  Now,  from  (2.11),  we  get  that 

B{e,e)  =—  f  e‘^{x,T)dx+  f  f  e^{x,t)dxdt— -  f  ey^{x,0)dx, 

2  Jq  Jo  Jo  ^  Jo 

and  since  eu{x,0)  =  tto,i(a:;)  -  «o.2(®)  and  B{e,e)  =  0,  by  the  error  equation,  we 

immediately  obtain  the  error  estimate  we  sought: 

-  el{x,T)dx+  /  el{x,t)dxdt^-  j  {uo,i{x)  -  uo,2{x)f  dx. 

2  Jq  Jo  Jo  ^  Jo 

To  prove  Theorem  2.2,  we  only  need  to  obtain  a  discrete  version  of  this  argument. 

The  discrete  case.  Since, 


Bh{'Wh,'Vh)  =0)  Vv/i(t)  eVhxVh,  Vt  e  (0,T), 
B/i(w,  v/i)  =  0,  Vv/i(t)  G  Vh  X  14,  VtG(0,T), 


by  (2.12a)  and  by  equations  (2.5a)  and  (2.5b),  respectively,  we  immediately  obtain 
our  error  equation: 

^h(e,V;i)=0,  Vv;i(t)  G  14  X  14,  VtG(0,T), 


where  e  =  w—Wh.  Now,  according  to  the  continuous  case  argument,  we  should  con¬ 
sider  next  the  quantity  B}i{b,  e);  however,  since  e  is  not  in  the  finite  element  space,  it 
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ismoreconvenienttoconsiderjB/i(P^(e),IP/i(e)),  where  IPh(e(t))  =  (P/i(ett(t)),  P/i(eg(t)) ) 
is  the  so-called  L^-projection  of  e(t)  into  the  finite  element  space  X  Vf.  The 
L^-projection  of  the  function  p  into  Vh,  f‘h{p),  is  defined  as  the  only  element  of  the 
finite  element  space  Vh  such  that 

/  {Fh{p){x)  -  p{x))  Vh{x)  dx  =  0,  'ivh^Vh-  (2.13) 

Jo 

Note  that,  in  fact  Uh{t  =  0)  =  Ph('i^o)>  by  (2.6c). 

Thus,  by  Lemma  2.3,  we  have 

Bh{Ph{e),Fh{e))=l  f  \Fh{eu{T)){x)\'^dx+  f  f  \Fh{eq{t)){x)\^dxdt 
^  Jo  Jo  Jo 

+  0r,c(lP/.(e)l)-i^‘  |P4e„(0))Wpir, 

and  since 

lPh(e«(0))  =  Ph(«o  -  uh{0))  =  Fh{uo)  -  Uh{0)  =  0, 
by  (2.6c)  and  (2.13),  and 

Bh(Fh{e),Fhie))  =  Bh{Fh{e)  -  e,Fhie))  =  ^/,(P,,(w)  -  w,Ph(e)), 

by  the  error  equation,  we  get 

1  T  1 

If  \Fh{eu{T)){x)\Ux+  f  f  \Fh{e,mx)\^dxdt 
^  Jo  Jo  Jo 

+  ©T,c([IPh(e)])  =  Bh{Fhi-w)-w,Fh{e)). 

(2.14) 

Note  that  since  in  our  continuous  model,  the  right-hand  side  is  zero,  we  expect  the 
term  ;S(P/i(w)  —  w,  P;i(e))  to  be  small. 

Estimating  the  right-hand  side.  To  show  that  this  is  so,  we  must  suitably 
treat  the  term  jB(Ph(w)  —  w,P;i(e)). 
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Lemma  2.4.  For  p  =  IP;i(w)  -  w,  we  have 


IPh(e))  <  2®t,c  (P)  +  ^^^  \Fhieqit)){x)  dx  dt 

+  {Axf^  r  C'a(t)  dt+{Ax)’^  r  C2{t)  (  /' 

Jo  Jo  I  Jo 


^hieu{t)){x)  {“^dx 


where 

Ci{t)  =  2  c|  I  ( Aa:|  ci2 1^  4 )  I  Ifc+i  +  a  (Aa;)^ |  «(t)  |, 

C2(i)  =  2cfc|  Vo|ci2  |tt(t)  |fc4.2  +  a(Aa;)^^^  ^M'“(0i|+2|- 

where  the  constants  Ck  and  dk  depend  solely  on  k,  and  k  —  k  except  when  the  grids 
are  uniform  and  k  is  even,  in  which  case  fc  =  fc  +  1. 


Note  how  Cii  appears  in  the  denominator  of  Ci{t).  However,  Ci{t)  remains 
bounded  as  cn  goes  to  zero  since  the  convective  numerical  flux  is  an  E-flux. 

To  prove  this  result,  we  will  need  the  following  auxiliary  lemmas.  We  denote  by 
I  'U'  Ik<*^+U(7)  integral  over  J  of  the  square  of  the  {k  +  l)-the  derivative  of  u. 

Lemma  2.5.  For  p  =  P;i(w)  —  w,  we  have 
|Puj+l/2l  ^  \ 

I  [Pu]j+l/2\<  Ck{Axf+'^/^\u\H(k+i)^j.^^^^), 

I  [Pq]j+i/2  I  <  CkVa{Ax)^+'^^^\u\jj(k+2)(j.^^^^), 

where  Jj+112  =  Ij  U  Ij+i,  the  constant  Ck  depends  solely  on  k,  and  k  =  k  except 
when  the  grids  are  uniform  and  k  is  even,  in  which  case  k  =  k  +  1. 


Proof.  The  two  last  inequalities  foUow  from  the  first  two  and  from  the  fact  that 
q  =  ■s/d  dxU.  The  two  first  inequalities  with  k  =  k  follow  from  the  definitions  of  pZ 
and  [Pu]  and  from  the  following  estimate: 

I  Ph('^)(®^l/2)  ~ '*^j+l/2  1  ^  2^kiAx)  ^  I  |7S'(*=+i)(J3.,.i/2)* 
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where  the  constant  Ck  depends  solely  on  k.  This  inequality  follows  from  the  fact 
that  “  %+i/2  =  0  when  n  is  a  polynomial  of  degree  k  and  from  a 

simple  application  of  the  Bramble-Hilbert  lemma;  see  Ciarlet  [11]. 

To  prove  the  inequalities  in  the  case  in  which  fc  =  fc  +  1,  we  only  need  to  show 
that  if  tt  is  a  polynomial  of  degree  fc  +  1  for  fc  even,  then  =  0.  It  is  clear  that  it 
is  enough  to  show  this  equality  for  the  particular  choice 

u{x)  =  {{x  -  Xj+i/2)/{Ax/2))’''^^ . 


To  prove  this,  we  recall  that  if  Pe  denotes  the  Legendre  polynomials  of  order  £: 
(i)  Pe{s)  Pm{s)  ds  =  (ii)  Pe{±l)  =  (±l)^  and  (iii)  Pe{s)  is  a  linear 

combination  of  odd  (even)  powers  of  s  for  odd  (even)  values  of  i.  Since  we  are 
assuming  that  the  grid  is  uniform,  Axj  =  Axj+i  =  Ax,  we  can  write,  by  (i). 


Ph{u){x)^  ^  ^^^^J^Pi{s)u{xj  +  ^Axs)ds'^Pe{ 


X  —  Xj 

Ax/2 


), 


for  X  E  Ij.  Hence,  for  our  particular  choice  of  u,  we  have 

K,+./2=5  E  ^  f  P,{s){(s-lf*'P,{l)  +  {s  +  lf*'Pt(-l)}dS 
^  o<e<k 

=  5  E  t  0  /*,  {(-l)'-"'"'  A(l)  +  P((-l)}  ds 


2  ^  ^  2 


=  l  E  ft +(-!)'}* 

^  0<£,i<k  ' 


by  (ii).  When  the  factor  {(-1)^+^“^  +  (-1)^}  is  different  from  zero,  \  k  +  l-i  +  £\ 
is  even  and  since  k  is  also  even,  j  i  f  |  is  odd.  In  this  case,  by  (iii). 


J  P^(s)s*ds  =  0, 


and  so  Puj+1/2  —  0-  This  completes  the  proof.  □ 

We  will  also  need  the  following  result  that  follows  from  a  simple  scaling  argument; 


see  Ciarlet  [11]. 
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Lemma  2.6.  We  have 


[¥h(p)]j+i/2  I  <  4  (Aa;)  ||  Fh{p)  ||i,2(a). 


where  4+1/2  =  4  4+1  constant  dk  depends  solely  on  k. 

We  are  now  ready  to  prove  Lemma  2.4. 


Proof  of  Lemma  2.4.  To  simplify  the  notation,  let  us  set  v/i  =  P/ie.  By  the  defini¬ 
tion  of  Bh{-,  ■),  we  have 


Bhijp^^h)  —  ['  f  4  4  4-  f  f  Pq(x,h)V}i^q{x^f)  dx  dt 

Jo  Jo  Jo  Jo 

rjn 

-  f  MP)i+l/2(4  [Vh(t)  4+1/2 

•'O  l<j<N 

-  f  X/  /  h(p(®)4)*4  Vh(a;,f)da;dt 

=  -  /  %)$-+l/2(4[vii(44+l/2Cfi, 


l<j<^ 


by  the  definition  of  the  L^ -projection  (2.13). 

Now,  recalling  that  p  =  {pu,PqY  and  that  v/i  =  {vu^VqY,  we  have 


h(p)*[v/i(t)]  =  {CPZ,-  Cll[Pu])[Vu] 

Jr{-Vap4-ci2[Pq])[vu] 
+  (-\/a^  +  Ci2  [Pu])  [Vq] 
=  01  -j-  02  p  ^3- 


By  Lemmas  2.5  and  2.6, 

I  ^1 1  <  Cfc  |«  lrrfc+i(j)  (I  c|  +  cn)  1  [wu]  |, 

I  02  I  ^  Cfc  dfc  (Aa?)^  (a  \  u  (A®)*  -t-  -x/o  |  C12  |  |  u  |ii'*:+2(j)))  ||  Vu  ||l2(j)) 

I  03  I  <  Ckdk  (AxY  (Va|'u|^fc+i(j)  (Aa;)^“^  +  1  C12  1 1 u |Hfc+i(jr)) )  \\vq  ||l2(j). 
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This  is  the  crucial  step  for  obtaining  our  error  estimates.  Note  that  the  treatment  of 
9 1  is  very  different  than  the  treatment  of  02  and  03.  The  reason  for  this  difference  is 
that  the  upper  bound  for  0i  can  be  controlled  by  the  form  ©T,c([vft])-  we  recall  that 
=  Fh{e).  This  is  not  the  case  for  the  upper  bound  for  02  because  ©t,c[v;i]  =  0 
if  c  =  0  nor  it  is  the  case  for  the  upper  bound  for  03  because  ©t,c[v;i]  does  not 
involve  the  jumps 

Thus,  after  a  suitable  application  of  Young’s  inequality  and  simple  algebraic 
manipulations,  we  get 

HpYiMt)]  <  \  cn  [Vu?  +  \\\V,  \  Ci(0  (Ax)^^  +  i  C2{t)  {Axf  II  Vu \W^jy 


Since 

J3h{p,Vh)<  [  Yj  |h(p)j+i/2(*)[vh(0b+i/2|dt, 

■^0  l<j<N 

and  since  J,+i/2  =  Ij  b)  Ij+u  the  result  follows  after  simple  applications  of  the 
Cauchy-Schwartz  inequality.  This  completes  the  proof.  □ 


Conclusion.  Combining  the  equation  (2.14)  with  the  estimate  of  Lemma  2.4, 
we  easily  obtain,  after  a  simple  application  of  Gronwall’s  lemma, 

^‘|Pk(e.(r))Wl^&  +  ^Y  I  {eq  (t) )  (x)  I  ^  dx  dt  +  ©t,c  ( (e)] )  | 


<  (Ax)'=  £  ^/^ dt  +  (Ax)''  £  C2{t)  I  I  Ph{eu{t)){x)  |2  dx  | 


1/2 


dt. 


Theorem  2.2  follows  easily  from  this  inequality.  Lemma  2.6,  and  from  the  following 
simple  approximation  result: 


\\p-Fh{p)  11^2(0, 1)  <  0fc(Ax)''+^  |p|^^(fe+i)(o,i) 
where  Qk  depends  solely  on  fc;  see  Ciarlet  [11]. 
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3.  The  LDG  methods  for  the  multi-dimensional  case.  In  this  section,  we 
consider  the  LDG  methods  for  the  following  convection-diffusion  model  problem 

dtu+  Y.  Y  o.iiin)d^.u)={)  in(0,T)x(0,l)^  (3.1a) 

\<i<d  l<j<d 

u{t  =  0)  =  «o,  on  (0, 1)'^,  (3.1b) 

with  periodic  boundary  conditions.  Essentially,  the  one-dimensional  case  and  the 
multidimensional  case  can  be  studied  in  exactly  the  same  way.  However,  there  are 
two  important  differences  that  deserve  explicit  discussion.  The  first  is  the  treatment 
of  the  matrix  of  entries  aij{u),  which  is  assumed  to  be  symmetric,  semipositive 
definite  and  the  introduction  of  the  variables  qi,  and  the  second  is  the  treatment  of 
arbitrary  meshes. 

To  define  the  LDG  method,  we  first  notice  that,  since  the  matrix  ay(tt)  is  as¬ 
sumed  to  be  symmetric  and  semipositive  definite,  there  exists  a  symmetric  matrix 
bij  (u)  such  that 

Y  hi{u)b£j{u).  (3.2) 

i<e<d 

Then  we  define  the  new  scalar  variables  qe  =  rewrite  the 

problem  (3.1)  as  follows: 

dtu+  Y  E  bie{u)qi)  =  0  in  (0, T)  x  (0, 1)'",  (3.3a) 

l<e<d 

qe-  Y  €  =  !,.. .d,  in  (0,T)  x  (0, 1)-",  (3.3b) 

l<j<d 

u{t  =  0)  =  uo,  on  (0, 1)^*,  (3.3c) 

where  gej{u)  =  bej{s)ds.  The  LDG  method  is  now  obtained  by  discretizing 
(3.3)  by  the  Discontinuous  Galerkin  method. 

We  follow  what  was  done  in  §2.  So,  we  set  w  =  {u,  q)*  =  (u,  9i,  •  ■  •  ,  qdY  and, 
for  each  *  =  1,  •  ■  •  ,d,  introduce  the  flux 

hi(w)  =  ( fi{u)  -  Y  qi,-gii{u),-  ■  •  ,  -gdi{u)  )^  (3.4) 

i<e<d 
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We  consider  triangulations  of  (0, 1)'^,  Ta®  =  {K},  made  of  non-overlapping  poly- 
hedra.  We  require  that  for  any  two  elements  /iT  and  K',  K  D  K  is  either  a  face 
e  of  both  K  and  K'  with  nonzero  {d  —  1)-Lebesgue  measure  |  e  |,  or  has  Hausdorff 
dimension  less  than  d  —  1.  We  denote  by  the  set  of  all  faces  e  of  the  border 
of  K  for  all  K  6  T^x-  The  diameter  of  K  is  denoted  by  Axk  and  the  maximum 
Axk,  for  K  E  T/S.X  is  denoted  by  Ax.  We  require,  for  the  sake  of  simplicity,  that 
the  triangulations  Tai;  be  regular,  that  is,  there  is  a  constant  independent  of  Ax 
such  that 

\/KeT^x, 

PK 

where  pK  denotes  the  diameter  of  the  maximum  ball  included  in  K. 

We  seek  an  approximation  Wh  =  {uh,  qhY  =  {^h,  Qhi,  ,  QhdY  to  w  such  that 
for  each  time  t  E  [0,T],  each  of  the  components  of  belong  to  the  finite  element 
space 

14  =  vjf  =  {t;  e  L\{0,  lY)  :  v\k  e  P\K)  yXE  r^xh  (3-5) 


where  P^{K)  denotes  the  space  of  polynomials  of  total  degree  at  most  k.  In  or¬ 
der  to  determine  the  approximate  solution  Wh,  we  proceed  exactly  as  in  the  one¬ 
dimensional  case.  This  time,  however,  the  integrals  are  made  on  each  element  K  of 
the  triangulation  T/^x-  We  obtain  the  following  weak  formulation  on  each  element 
K  of  the  triangulation  Ta*: 


/  dtUh{x,t)vh,u{x)dx-  Y]  /  hiu{-Wh{x,t))dc,^Vh,u{x)dx 


+  f  K{wh,ndK){x,t)vh,u{x)dr{x)  =0,  yvh,u^P^{K), 
JdK 


For  £  =  1,  -  ••  ,d  : 


Vh,qe{x)dx 


{Wh{x,  t))  Vh,qi{^)  dx 

i<i<d 


+  f  hge{'^h,ndK){x,t)Vh,qt{x)dT{x)  =0,  ^  Vk,„  e  P^K), 

JdK 


(3.6a) 


(3.6b) 
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/  uh{x,0)vh,i{x)dx  =  /  uo{x) Vh,i{x)dx,  W Vh,i  E  P’'{K), 

Jk  Jk  (3.6c) 

where  xiqk  denotes  the  outward  unit  normal  to  the  element  K  at  x  E  dK.  It 
remains  to  choose  the  numerical  flux  (hu,  ,  Kd)*  =  h  =  h(w^,  nai;c)(a;,  t). 

As  in  the  one-dimensional  case,  we  require  that  the  fluxes  h  be  of  the  form 

h{-Wh,neK)ix)  =  na^r), 

where  is  the  limit  at  x  taken  from  the  interior  of  K  and  Wh{x^^^^)  the 

limit  at  X  from  the  exterior  of  K,  and  consider  fluxes  that  (i)  are  locally  Lipschitz, 
conservative,  that  is, 


and  consistent  with  the  flux 

^  ^  hj  ndK,if 
l<i<d 

(ii)  allow  for  a  local  resolution  of  each  component  of  in  terms  of  Uh  only,  (iii) 
reduce  to  an  E-flux  when  a(-)  =  0,  and  that  (iv)  enforce  the  L^-stability  of  the 
method. 

Again,  we  write  our  numerical  flux  as  the  sum  of  a  convective  flux  and  a  diffusive 
flux: 


h  —  hconu  “b  ^diffy 


where  the  convective  flux  is  given  by 

hcont;(w“,w+;n)  =  (/(«",  ^t+;n),0)^ 

where  /(u“,u+;  n)  is  any  locally  Lipschitz  E-flux  which  is  conservative  and  consis¬ 
tent  with  the  nonlinearity 

fi{u)ni, 

i<i<d 
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and  the  diffusive  flux  hrfi//(w  ,  w+;  n)  is  given  by 

(-  ^  ^  Pid('u)ni  )*-Cdi// [w], 

l<i,e<d  l<i<<i  l<i<d 

where 


( ° 
-Cl2 


^diff  = 


-Cl3 


Cl2  Ci3  •  •  •  Cirf  \ 

0  0  0 

0  0  ■■■  0 


V  —cid  0  0  ■  •  ■  0  / 


cij  =  cij(w  jW"^)  is  locally  Lipschitz  for  j  =  1,  ■  •  ■  ,  d, 


cij  =  0  when  a(-)  =  0  for  j  =  1,  ■  -  ■  ,  d. 


We  claim  that  this  flux  satisfies  the  properties  (i)  to  (iv). 

To  prove  that  properties  (i)  to  (iii)  are  satisfied  is  now  a  simple  exercise.  To  see 
that  the  property  (iv)  is  satisfied,  we  first  rewrite  the  flux  h  in  the  following  way: 

l<i,i<d  l<i<d  l<i<d 


where 


/  Cii 

I  -Cl2 


c  = 


-C13 


Cl2  Ci3 
0  0 

0  0 


Cld\ 

0 


\—cid  0  0 


0  / 

/(«",tt+;n)  Y 


where  <f>i{u)  =  fi{s)  ds.  Since  /(-,  •;n)  is  an  E-flux, 


Cll 


(  S  ,ti'^;n))ds  >0, 

l<i<d 


and  so  the  matrix  C  is  semipositive  definite.  The  property  (iv)  follows  from  this 
fact  and  from  the  following  multidimensional  version  of  Proposition  2.1. 
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Proposition  3.1.  (L^-stability)  We  have, 


where 


0T,c([wh])=  /  V  [  [wh{x,t)YC[wh{x,t)]dT{x)dt. 

We  can  also  prove  the  following  error  estimate.  We  denote  the  integral  over 
(0, 1)*^  of  the  sum  of  the  squares  of  all  the  derivatives  of  order  {k  +  1)  of  «  by 


\u 


fc+i- 


Theorem  3.2.  (L^-error  estimate)  Let  e  he  the  approximation  error  w—Wh.  Then 
we  have,  for  arbitrary,  regular  grids, 

I  [  \eu{x,T)\'^dx+  [  [  leg(a;,t)  pda:dt  +  ©T,c([e])|  <C{Ax)^, 

I  J(o,iy  Jo  J{o,iy  J 


where  (7  =  (7(A:,  | « |fc+i,  |  n  |fc+2)-  In  the  purely  hyperbolic  case  aij  =  0,  the  constant 
C  is  of  order  In  the  purely  parabolic  case  c  =  0,  the  constant  C  is  of  order 

Ax  for  even  values  of  k  and  of  order  1  otherwise  for  Cartesian  products  of  uniform 
grids  and  for  C  identically  zero  provided  that  the  local  spaces  are  used  instead 
of  the  spaces  ,  where  is  the  space  of  tensor  products  of  one  dimensional 
polynomials  of  degree  k. 


The  algebraic  manipulations  needed  to  prove  this  result  are  a  straightforward 
extension  to  the  multidimensional  case  of  the  manipulations  of  the  proof  of  the 
corresponding  one-dimensional  result,  Theorem  2.2.  The  approximation  properties 
of  the  finite  element  spaces  Vh  that  extend  the  results  of  Lemmas  2.5  and  2.6  are 
the  following.  Let  e  denote  a  face  of  the  element  K  and  let  us  denote  by  the 
element  sharing  the  face  e  with  K,  then 

-  wllL2(e)  <  -CkiAx)’^'^^^'^  \ 
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where  Vh{u)'^  is  either  the  value  of  Fh{u)  at  e  from  the  interior  of  K  or  from  its 
exterior,  and 

II  [P/i(p)]  WL^ie)  <  dk  ||  Fh{p)  IU2(jcuii'e)) 

where  [Ph(p)]  denotes  the  jump  at  e  of  Ph(p)-  Finally,  we  also  use  the  following 
result: 


||p- P/i(p)  IU2(o,i)<i  <  |p|;f{fc+i)(0,i)d, 

All  these  approximation  results  can  be  found  in  Ciarlet  [11],  for  example. 

4.  Concluding  remarks.  In  this  paper,  we  have  considered  the  so-called  LDG 
methods  for  convection-diffusion  problems.  For  scalar  problems  in  multidimensions, 
we  have  shown  that  they  are  L^-stable  and  that  in  the  linear  case,  they  are  of  order 
k  if  polynomials  of  order  k  are  used.  We  have  also  shown  that  this  estimate  is  sharp 
and  have  displayed  the  strong  dependence  of  the  order  of  convergence  of  the  LDG 
methods  on  the  choice  of  the  numerical  fluxes. 

The  LDG  methods  for  multidimensional  systems,  like  for  example  the  compress¬ 
ible  Navier-Stokes  equations  and  the  equations  of  the  hydrodynamic  model  for 
semiconductor  device  simulation,  can  be  easily  defined  by  simply  applying  the  pro¬ 
cedure  described  for  the  multidimensional  scalar  case  to  each  component  of  u.  In 
practice,  especially  for  viscous  terms  which  are  not  symmetric  but  still  semipos¬ 
itive  definite,  such  as  for  the  compressible  Navier-Stokes  equations,  we  can  use 
q  =  (dx^  u,  ...,dx^u)  as  the  auxiliary  variables.  Although  with  this  choice,  the  bi¬ 
stability  result  will  not  be  available  theoretically,  this  would  not  cause  any  problem 
in  practical  implementation;  see  Bassi  and  Rebay  [5]  and  Bassi  et  al  [6] . 

The  main  advantage  of  these  methods  is  their  extremely  high  parallelizabil- 
ity  and  their  high-order  accuracy  which  render  them  suitable  for  computations 
of  convection-dominated  flows.  Indeed,  although  the  LDG  method  have  a  large 
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amount  of  degrees  of  freedom  per  element,  and  hence  more  computations  per  ele¬ 
ment  are  necessary,  its  extremely  local  domain  of  dependency  allows  a  very  efficient 
parallelization  that  by  far  compensates  for  the  extra  amount  of  local  computations. 
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